home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zunk2.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  27.1 KB  |  647 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((zeror 0.0)
  12.       (zeroi 0.0)
  13.       (coner 1.0)
  14.       (cr1r 1.0)
  15.       (cr1i 1.7320508075688772)
  16.       (cr2r -0.5)
  17.       (cr2i -0.8660254037844388)
  18.       (hpi 1.5707963267948966)
  19.       (pi_ 3.141592653589793)
  20.       (aic 1.2655121234846454)
  21.       (cipr (make-array 4 :element-type 'double-float))
  22.       (cipi (make-array 4 :element-type 'double-float)))
  23.   (declare (type (simple-array double-float (4)) cipi cipr)
  24.            (type double-float aic pi_ hpi cr2i cr2r cr1i cr1r coner zeroi
  25.             zeror))
  26.   (f2cl-lib:fset (f2cl-lib:fref cipr (1) ((1 4))) 1.0)
  27.   (f2cl-lib:fset (f2cl-lib:fref cipi (1) ((1 4))) 0.0)
  28.   (f2cl-lib:fset (f2cl-lib:fref cipr (2) ((1 4))) 0.0)
  29.   (f2cl-lib:fset (f2cl-lib:fref cipi (2) ((1 4))) -1.0)
  30.   (f2cl-lib:fset (f2cl-lib:fref cipr (3) ((1 4))) -1.0)
  31.   (f2cl-lib:fset (f2cl-lib:fref cipi (3) ((1 4))) 0.0)
  32.   (f2cl-lib:fset (f2cl-lib:fref cipr (4) ((1 4))) 0.0)
  33.   (f2cl-lib:fset (f2cl-lib:fref cipi (4) ((1 4))) 1.0)
  34.   (defun zunk2 (zr zi fnu kode mr n yr yi nz tol elim alim)
  35.     (declare (type (simple-array double-float (*)) yr yi)
  36.              (type f2cl-lib:integer4 kode mr n nz)
  37.              (type double-float zr zi fnu tol elim alim))
  38.     (prog ((bry (make-array 3 :element-type 'double-float))
  39.            (asumr (make-array 2 :element-type 'double-float))
  40.            (asumi (make-array 2 :element-type 'double-float))
  41.            (bsumr (make-array 2 :element-type 'double-float))
  42.            (bsumi (make-array 2 :element-type 'double-float))
  43.            (phir (make-array 2 :element-type 'double-float))
  44.            (phii (make-array 2 :element-type 'double-float))
  45.            (argr (make-array 2 :element-type 'double-float))
  46.            (argi (make-array 2 :element-type 'double-float))
  47.            (zeta1r (make-array 2 :element-type 'double-float))
  48.            (zeta1i (make-array 2 :element-type 'double-float))
  49.            (zeta2r (make-array 2 :element-type 'double-float))
  50.            (zeta2i (make-array 2 :element-type 'double-float))
  51.            (cyr (make-array 2 :element-type 'double-float))
  52.            (cyi (make-array 2 :element-type 'double-float))
  53.            (cssr (make-array 3 :element-type 'double-float))
  54.            (csrr (make-array 3 :element-type 'double-float)) (i 0) (ib 0)
  55.            (iflag 0) (ifn 0) (il 0) (in 0) (inu 0) (iuf 0) (k 0) (kdflg 0)
  56.            (kflag 0) (kk 0) (nai 0) (ndai 0) (nw 0) (idum 0) (j 0) (ipard 0)
  57.            (ic 0) (aarg 0.0) (aii 0.0) (air 0.0) (ang 0.0) (aphi 0.0)
  58.            (argdi 0.0) (argdr 0.0) (asc 0.0) (ascle 0.0) (asumdi 0.0)
  59.            (asumdr 0.0) (bsumdi 0.0) (bsumdr 0.0) (car 0.0) (cki 0.0) (ckr 0.0)
  60.            (crsc 0.0) (cscl 0.0) (csgni 0.0) (csi 0.0) (cspni 0.0) (cspnr 0.0)
  61.            (csr 0.0) (c1i 0.0) (c1r 0.0) (c2i 0.0) (c2m 0.0) (c2r 0.0)
  62.            (daii 0.0) (dair 0.0) (fmr 0.0) (fn 0.0) (fnf 0.0) (phidi 0.0)
  63.            (phidr 0.0) (pti 0.0) (ptr 0.0) (rast 0.0) (razr 0.0) (rs1 0.0)
  64.            (rzi 0.0) (rzr 0.0) (sar 0.0) (sgn 0.0) (sti 0.0) (str 0.0)
  65.            (s1i 0.0) (s1r 0.0) (s2i 0.0) (s2r 0.0) (yy 0.0) (zbi 0.0) (zbr 0.0)
  66.            (zet1di 0.0) (zet1dr 0.0) (zet2di 0.0) (zet2dr 0.0) (zni 0.0)
  67.            (znr 0.0) (zri 0.0) (zrr 0.0))
  68.       (declare (type (simple-array double-float (3)) cssr csrr bry)
  69.                (type (simple-array double-float (2)) zeta2r zeta2i zeta1r
  70.                 zeta1i phir phii cyr cyi bsumr bsumi asumr asumi argr argi)
  71.                (type double-float zrr zri znr zni zet2dr zet2di zet1dr zet1di
  72.                 zbr zbi yy s2r s2i s1r s1i str sti sgn sar rzr rzi rs1 razr
  73.                 rast ptr pti phidr phidi fnf fn fmr dair daii c2r c2m c2i c1r
  74.                 c1i csr cspnr cspni csi csgni cscl crsc ckr cki car bsumdr
  75.                 bsumdi asumdr asumdi ascle asc argdr argdi aphi ang air aii
  76.                 aarg)
  77.                (type f2cl-lib:integer4 ic ipard j idum nw ndai nai kk kflag
  78.                 kdflg k iuf inu in il ifn iflag ib i))
  79.       (setf kdflg 1)
  80.       (setf nz 0)
  81.       (setf cscl (/ 1.0 tol))
  82.       (setf crsc tol)
  83.       (f2cl-lib:fset (f2cl-lib:fref cssr (1) ((1 3))) cscl)
  84.       (f2cl-lib:fset (f2cl-lib:fref cssr (2) ((1 3))) coner)
  85.       (f2cl-lib:fset (f2cl-lib:fref cssr (3) ((1 3))) crsc)
  86.       (f2cl-lib:fset (f2cl-lib:fref csrr (1) ((1 3))) crsc)
  87.       (f2cl-lib:fset (f2cl-lib:fref csrr (2) ((1 3))) coner)
  88.       (f2cl-lib:fset (f2cl-lib:fref csrr (3) ((1 3))) cscl)
  89.       (f2cl-lib:fset (f2cl-lib:fref bry (1) ((1 3)))
  90.                      (/ (* 1000.0 (f2cl-lib:d1mach 1)) tol))
  91.       (f2cl-lib:fset (f2cl-lib:fref bry (2) ((1 3)))
  92.                      (/ 1.0 (f2cl-lib:fref bry (1) ((1 3)))))
  93.       (f2cl-lib:fset (f2cl-lib:fref bry (3) ((1 3))) (f2cl-lib:d1mach 2))
  94.       (setf zrr zr)
  95.       (setf zri zi)
  96.       (if (>= zr 0.0) (go label10))
  97.       (setf zrr (- zr))
  98.       (setf zri (- zi))
  99.      label10
  100.       (setf yy zri)
  101.       (setf znr zri)
  102.       (setf zni (- zrr))
  103.       (setf zbr zrr)
  104.       (setf zbi zri)
  105.       (setf inu (f2cl-lib:int fnu))
  106.       (setf fnf (- fnu inu))
  107.       (setf ang (* (- hpi) fnf))
  108.       (setf car (cos ang))
  109.       (setf sar (sin ang))
  110.       (setf c2r (* hpi sar))
  111.       (setf c2i (* (- hpi) car))
  112.       (setf kk (f2cl-lib:int-add (mod inu 4) 1))
  113.       (setf str
  114.               (- (* c2r (f2cl-lib:fref cipr (kk) ((1 4))))
  115.                  (* c2i (f2cl-lib:fref cipi (kk) ((1 4))))))
  116.       (setf sti
  117.               (+ (* c2r (f2cl-lib:fref cipi (kk) ((1 4))))
  118.                  (* c2i (f2cl-lib:fref cipr (kk) ((1 4))))))
  119.       (setf csr (- (* cr1r str) (* cr1i sti)))
  120.       (setf csi (+ (* cr1r sti) (* cr1i str)))
  121.       (if (> yy 0.0) (go label20))
  122.       (setf znr (- znr))
  123.       (setf zbi (- zbi))
  124.      label20
  125.       (setf j 2)
  126.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  127.                     ((> i n) nil)
  128.         (tagbody
  129.           (setf j (f2cl-lib:int-sub 3 j))
  130.           (setf fn (+ fnu (f2cl-lib:int-sub i 1)))
  131.           (multiple-value-bind
  132.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
  133.                var-10 var-11 var-12 var-13 var-14 var-15 var-16)
  134.               (zunhj znr zni fn 0 tol (f2cl-lib:fref phir (j) ((1 2)))
  135.                (f2cl-lib:fref phii (j) ((1 2)))
  136.                (f2cl-lib:fref argr (j) ((1 2)))
  137.                (f2cl-lib:fref argi (j) ((1 2)))
  138.                (f2cl-lib:fref zeta1r (j) ((1 2)))
  139.                (f2cl-lib:fref zeta1i (j) ((1 2)))
  140.                (f2cl-lib:fref zeta2r (j) ((1 2)))
  141.                (f2cl-lib:fref zeta2i (j) ((1 2)))
  142.                (f2cl-lib:fref asumr (j) ((1 2)))
  143.                (f2cl-lib:fref asumi (j) ((1 2)))
  144.                (f2cl-lib:fref bsumr (j) ((1 2)))
  145.                (f2cl-lib:fref bsumi (j) ((1 2))))
  146.             (declare (ignore var-0 var-1 var-2 var-3 var-4))
  147.             (f2cl-lib:fset (f2cl-lib:fref phir (j) ((1 2))) var-5)
  148.             (f2cl-lib:fset (f2cl-lib:fref phii (j) ((1 2))) var-6)
  149.             (f2cl-lib:fset (f2cl-lib:fref argr (j) ((1 2))) var-7)
  150.             (f2cl-lib:fset (f2cl-lib:fref argi (j) ((1 2))) var-8)
  151.             (f2cl-lib:fset (f2cl-lib:fref zeta1r (j) ((1 2))) var-9)
  152.             (f2cl-lib:fset (f2cl-lib:fref zeta1i (j) ((1 2))) var-10)
  153.             (f2cl-lib:fset (f2cl-lib:fref zeta2r (j) ((1 2))) var-11)
  154.             (f2cl-lib:fset (f2cl-lib:fref zeta2i (j) ((1 2))) var-12)
  155.             (f2cl-lib:fset (f2cl-lib:fref asumr (j) ((1 2))) var-13)
  156.             (f2cl-lib:fset (f2cl-lib:fref asumi (j) ((1 2))) var-14)
  157.             (f2cl-lib:fset (f2cl-lib:fref bsumr (j) ((1 2))) var-15)
  158.             (f2cl-lib:fset (f2cl-lib:fref bsumi (j) ((1 2))) var-16))
  159.           (if (= kode 1) (go label30))
  160.           (setf str (+ zbr (f2cl-lib:fref zeta2r (j) ((1 2)))))
  161.           (setf sti (+ zbi (f2cl-lib:fref zeta2i (j) ((1 2)))))
  162.           (setf rast (/ fn (zabs str sti)))
  163.           (setf str (* str rast rast))
  164.           (setf sti (* (- sti) rast rast))
  165.           (setf s1r (- (f2cl-lib:fref zeta1r (j) ((1 2))) str))
  166.           (setf s1i (- (f2cl-lib:fref zeta1i (j) ((1 2))) sti))
  167.           (go label40)
  168.          label30
  169.           (setf s1r
  170.                   (- (f2cl-lib:fref zeta1r (j) ((1 2)))
  171.                      (f2cl-lib:fref zeta2r (j) ((1 2)))))
  172.           (setf s1i
  173.                   (- (f2cl-lib:fref zeta1i (j) ((1 2)))
  174.                      (f2cl-lib:fref zeta2i (j) ((1 2)))))
  175.          label40
  176.           (setf rs1 s1r)
  177.           (if (> (abs rs1) elim) (go label70))
  178.           (if (= kdflg 1) (setf kflag 2))
  179.           (if (< (abs rs1) alim) (go label50))
  180.           (setf aphi
  181.                   (zabs (f2cl-lib:fref phir (j) ((1 2)))
  182.                    (f2cl-lib:fref phii (j) ((1 2)))))
  183.           (setf aarg
  184.                   (zabs (f2cl-lib:fref argr (j) ((1 2)))
  185.                    (f2cl-lib:fref argi (j) ((1 2)))))
  186.           (setf rs1
  187.                   (- (+ rs1 (f2cl-lib:flog aphi))
  188.                      (* 0.25 (f2cl-lib:flog aarg))
  189.                      aic))
  190.           (if (> (abs rs1) elim) (go label70))
  191.           (if (= kdflg 1) (setf kflag 1))
  192.           (if (< rs1 0.0) (go label50))
  193.           (if (= kdflg 1) (setf kflag 3))
  194.          label50
  195.           (setf c2r
  196.                   (- (* (f2cl-lib:fref argr (j) ((1 2))) cr2r)
  197.                      (* (f2cl-lib:fref argi (j) ((1 2))) cr2i)))
  198.           (setf c2i
  199.                   (+ (* (f2cl-lib:fref argr (j) ((1 2))) cr2i)
  200.                      (* (f2cl-lib:fref argi (j) ((1 2))) cr2r)))
  201.           (multiple-value-bind
  202.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  203.               (zairy c2r c2i 0 2 air aii nai idum)
  204.             (declare (ignore var-0 var-1 var-2 var-3))
  205.             (setf air var-4)
  206.             (setf aii var-5)
  207.             (setf nai var-6)
  208.             (setf idum var-7))
  209.           (multiple-value-bind
  210.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  211.               (zairy c2r c2i 1 2 dair daii ndai idum)
  212.             (declare (ignore var-0 var-1 var-2 var-3))
  213.             (setf dair var-4)
  214.             (setf daii var-5)
  215.             (setf ndai var-6)
  216.             (setf idum var-7))
  217.           (setf str
  218.                   (- (* dair (f2cl-lib:fref bsumr (j) ((1 2))))
  219.                      (* daii (f2cl-lib:fref bsumi (j) ((1 2))))))
  220.           (setf sti
  221.                   (+ (* dair (f2cl-lib:fref bsumi (j) ((1 2))))
  222.                      (* daii (f2cl-lib:fref bsumr (j) ((1 2))))))
  223.           (setf ptr (- (* str cr2r) (* sti cr2i)))
  224.           (setf pti (+ (* str cr2i) (* sti cr2r)))
  225.           (setf str
  226.                   (+ ptr
  227.                      (- (* air (f2cl-lib:fref asumr (j) ((1 2))))
  228.                         (* aii (f2cl-lib:fref asumi (j) ((1 2)))))))
  229.           (setf sti
  230.                   (+ pti
  231.                      (+ (* air (f2cl-lib:fref asumi (j) ((1 2))))
  232.                         (* aii (f2cl-lib:fref asumr (j) ((1 2)))))))
  233.           (setf ptr
  234.                   (- (* str (f2cl-lib:fref phir (j) ((1 2))))
  235.                      (* sti (f2cl-lib:fref phii (j) ((1 2))))))
  236.           (setf pti
  237.                   (+ (* str (f2cl-lib:fref phii (j) ((1 2))))
  238.                      (* sti (f2cl-lib:fref phir (j) ((1 2))))))
  239.           (setf s2r (- (* ptr csr) (* pti csi)))
  240.           (setf s2i (+ (* ptr csi) (* pti csr)))
  241.           (setf str (* (exp s1r) (f2cl-lib:fref cssr (kflag) ((1 3)))))
  242.           (setf s1r (* str (cos s1i)))
  243.           (setf s1i (* str (sin s1i)))
  244.           (setf str (- (* s2r s1r) (* s2i s1i)))
  245.           (setf s2i (+ (* s1r s2i) (* s2r s1i)))
  246.           (setf s2r str)
  247.           (if (/= kflag 1) (go label60))
  248.           (multiple-value-bind
  249.               (var-0 var-1 var-2 var-3 var-4)
  250.               (zuchk s2r s2i nw (f2cl-lib:fref bry (1) ((1 3))) tol)
  251.             (declare (ignore var-0 var-1 var-3 var-4))
  252.             (setf nw var-2))
  253.           (if (/= nw 0) (go label70))
  254.          label60
  255.           (if (<= yy 0.0) (setf s2i (- s2i)))
  256.           (f2cl-lib:fset (f2cl-lib:fref cyr (kdflg) ((1 2))) s2r)
  257.           (f2cl-lib:fset (f2cl-lib:fref cyi (kdflg) ((1 2))) s2i)
  258.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n)))
  259.                          (* s2r (f2cl-lib:fref csrr (kflag) ((1 3)))))
  260.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n)))
  261.                          (* s2i (f2cl-lib:fref csrr (kflag) ((1 3)))))
  262.           (setf str csi)
  263.           (setf csi (- csr))
  264.           (setf csr str)
  265.           (if (= kdflg 2) (go label85))
  266.           (setf kdflg 2)
  267.           (go label80)
  268.          label70
  269.           (if (> rs1 0.0) (go label320))
  270.           (if (< zr 0.0) (go label320))
  271.           (setf kdflg 1)
  272.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) zeror)
  273.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) zeroi)
  274.           (setf nz (f2cl-lib:int-add nz 1))
  275.           (setf str csi)
  276.           (setf csi (- csr))
  277.           (setf csr str)
  278.           (if (= i 1) (go label80))
  279.           (if
  280.            (and (= (f2cl-lib:fref yr ((f2cl-lib:int-sub i 1)) ((1 n))) zeror)
  281.                 (= (f2cl-lib:fref yi ((f2cl-lib:int-sub i 1)) ((1 n))) zeroi))
  282.            (go label80))
  283.           (f2cl-lib:fset (f2cl-lib:fref yr ((f2cl-lib:int-sub i 1)) ((1 n)))
  284.                          zeror)
  285.           (f2cl-lib:fset (f2cl-lib:fref yi ((f2cl-lib:int-sub i 1)) ((1 n)))
  286.                          zeroi)
  287.           (setf nz (f2cl-lib:int-add nz 1))
  288.          label80))
  289.       (setf i n)
  290.      label85
  291.       (setf razr (/ 1.0 (zabs zrr zri)))
  292.       (setf str (* zrr razr))
  293.       (setf sti (* (- zri) razr))
  294.       (setf rzr (* (+ str str) razr))
  295.       (setf rzi (* (+ sti sti) razr))
  296.       (setf ckr (* fn rzr))
  297.       (setf cki (* fn rzi))
  298.       (setf ib (f2cl-lib:int-add i 1))
  299.       (if (< n ib) (go label180))
  300.       (setf fn (+ fnu (f2cl-lib:int-sub n 1)))
  301.       (setf ipard 1)
  302.       (if (/= mr 0) (setf ipard 0))
  303.       (multiple-value-bind
  304.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  305.            var-11 var-12 var-13 var-14 var-15 var-16)
  306.           (zunhj znr zni fn ipard tol phidr phidi argdr argdi zet1dr zet1di
  307.            zet2dr zet2di asumdr asumdi bsumdr bsumdi)
  308.         (declare (ignore var-0 var-1 var-2 var-3 var-4))
  309.         (setf phidr var-5)
  310.         (setf phidi var-6)
  311.         (setf argdr var-7)
  312.         (setf argdi var-8)
  313.         (setf zet1dr var-9)
  314.         (setf zet1di var-10)
  315.         (setf zet2dr var-11)
  316.         (setf zet2di var-12)
  317.         (setf asumdr var-13)
  318.         (setf asumdi var-14)
  319.         (setf bsumdr var-15)
  320.         (setf bsumdi var-16))
  321.       (if (= kode 1) (go label90))
  322.       (setf str (+ zbr zet2dr))
  323.       (setf sti (+ zbi zet2di))
  324.       (setf rast (/ fn (zabs str sti)))
  325.       (setf str (* str rast rast))
  326.       (setf sti (* (- sti) rast rast))
  327.       (setf s1r (- zet1dr str))
  328.       (setf s1i (- zet1di sti))
  329.       (go label100)
  330.      label90
  331.       (setf s1r (- zet1dr zet2dr))
  332.       (setf s1i (- zet1di zet2di))
  333.      label100
  334.       (setf rs1 s1r)
  335.       (if (> (abs rs1) elim) (go label105))
  336.       (if (< (abs rs1) alim) (go label120))
  337.       (setf aphi (zabs phidr phidi))
  338.       (setf rs1 (+ rs1 (f2cl-lib:flog aphi)))
  339.       (if (< (abs rs1) elim) (go label120))
  340.      label105
  341.       (if (> rs1 0.0) (go label320))
  342.       (if (< zr 0.0) (go label320))
  343.       (setf nz n)
  344.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  345.                     ((> i n) nil)
  346.         (tagbody
  347.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) zeror)
  348.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) zeroi)
  349.          label106))
  350.       (go end_label)
  351.      label120
  352.       (setf s1r (f2cl-lib:fref cyr (1) ((1 2))))
  353.       (setf s1i (f2cl-lib:fref cyi (1) ((1 2))))
  354.       (setf s2r (f2cl-lib:fref cyr (2) ((1 2))))
  355.       (setf s2i (f2cl-lib:fref cyi (2) ((1 2))))
  356.       (setf c1r (f2cl-lib:fref csrr (kflag) ((1 3))))
  357.       (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
  358.       (f2cl-lib:fdo (i ib (f2cl-lib:int-add i 1))
  359.                     ((> i n) nil)
  360.         (tagbody
  361.           (setf c2r s2r)
  362.           (setf c2i s2i)
  363.           (setf s2r (+ (- (* ckr c2r) (* cki c2i)) s1r))
  364.           (setf s2i (+ (* ckr c2i) (* cki c2r) s1i))
  365.           (setf s1r c2r)
  366.           (setf s1i c2i)
  367.           (setf ckr (+ ckr rzr))
  368.           (setf cki (+ cki rzi))
  369.           (setf c2r (* s2r c1r))
  370.           (setf c2i (* s2i c1r))
  371.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) c2r)
  372.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) c2i)
  373.           (if (>= kflag 3) (go label130))
  374.           (setf str (coerce (abs c2r) 'double-float))
  375.           (setf sti (coerce (abs c2i) 'double-float))
  376.           (setf c2m (max str sti))
  377.           (if (<= c2m ascle) (go label130))
  378.           (setf kflag (f2cl-lib:int-add kflag 1))
  379.           (setf ascle (f2cl-lib:fref bry (kflag) ((1 3))))
  380.           (setf s1r (* s1r c1r))
  381.           (setf s1i (* s1i c1r))
  382.           (setf s2r c2r)
  383.           (setf s2i c2i)
  384.           (setf s1r (* s1r (f2cl-lib:fref cssr (kflag) ((1 3)))))
  385.           (setf s1i (* s1i (f2cl-lib:fref cssr (kflag) ((1 3)))))
  386.           (setf s2r (* s2r (f2cl-lib:fref cssr (kflag) ((1 3)))))
  387.           (setf s2i (* s2i (f2cl-lib:fref cssr (kflag) ((1 3)))))
  388.           (setf c1r (f2cl-lib:fref csrr (kflag) ((1 3))))
  389.          label130))
  390.      label180
  391.       (if (= mr 0) (go end_label))
  392.       (setf nz 0)
  393.       (setf fmr (coerce (the f2cl-lib:integer4 mr) 'double-float))
  394.       (setf sgn (coerce (- (f2cl-lib:dsign pi_ fmr)) 'double-float))
  395.       (setf csgni sgn)
  396.       (if (<= yy 0.0) (setf csgni (- csgni)))
  397.       (setf ifn (f2cl-lib:int-sub (f2cl-lib:int-add inu n) 1))
  398.       (setf ang (* fnf sgn))
  399.       (setf cspnr (cos ang))
  400.       (setf cspni (sin ang))
  401.       (if (= (mod ifn 2) 0) (go label190))
  402.       (setf cspnr (- cspnr))
  403.       (setf cspni (- cspni))
  404.      label190
  405.       (setf csr (* sar csgni))
  406.       (setf csi (* car csgni))
  407.       (setf in (f2cl-lib:int-add (mod ifn 4) 1))
  408.       (setf c2r (f2cl-lib:fref cipr (in) ((1 4))))
  409.       (setf c2i (f2cl-lib:fref cipi (in) ((1 4))))
  410.       (setf str (+ (* csr c2r) (* csi c2i)))
  411.       (setf csi (+ (* (- csr) c2i) (* csi c2r)))
  412.       (setf csr str)
  413.       (setf asc (f2cl-lib:fref bry (1) ((1 3))))
  414.       (setf iuf 0)
  415.       (setf kk n)
  416.       (setf kdflg 1)
  417.       (setf ib (f2cl-lib:int-sub ib 1))
  418.       (setf ic (f2cl-lib:int-sub ib 1))
  419.       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  420.                     ((> k n) nil)
  421.         (tagbody
  422.           (setf fn (+ fnu (f2cl-lib:int-sub kk 1)))
  423.           (if (> n 2) (go label175))
  424.          label172
  425.           (setf phidr (f2cl-lib:fref phir (j) ((1 2))))
  426.           (setf phidi (f2cl-lib:fref phii (j) ((1 2))))
  427.           (setf argdr (f2cl-lib:fref argr (j) ((1 2))))
  428.           (setf argdi (f2cl-lib:fref argi (j) ((1 2))))
  429.           (setf zet1dr (f2cl-lib:fref zeta1r (j) ((1 2))))
  430.           (setf zet1di (f2cl-lib:fref zeta1i (j) ((1 2))))
  431.           (setf zet2dr (f2cl-lib:fref zeta2r (j) ((1 2))))
  432.           (setf zet2di (f2cl-lib:fref zeta2i (j) ((1 2))))
  433.           (setf asumdr (f2cl-lib:fref asumr (j) ((1 2))))
  434.           (setf asumdi (f2cl-lib:fref asumi (j) ((1 2))))
  435.           (setf bsumdr (f2cl-lib:fref bsumr (j) ((1 2))))
  436.           (setf bsumdi (f2cl-lib:fref bsumi (j) ((1 2))))
  437.           (setf j (f2cl-lib:int-sub 3 j))
  438.           (go label210)
  439.          label175
  440.           (if (and (= kk n) (< ib n)) (go label210))
  441.           (if (or (= kk ib) (= kk ic)) (go label172))
  442.           (multiple-value-bind
  443.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
  444.                var-10 var-11 var-12 var-13 var-14 var-15 var-16)
  445.               (zunhj znr zni fn 0 tol phidr phidi argdr argdi zet1dr zet1di
  446.                zet2dr zet2di asumdr asumdi bsumdr bsumdi)
  447.             (declare (ignore var-0 var-1 var-2 var-3 var-4))
  448.             (setf phidr var-5)
  449.             (setf phidi var-6)
  450.             (setf argdr var-7)
  451.             (setf argdi var-8)
  452.             (setf zet1dr var-9)
  453.             (setf zet1di var-10)
  454.             (setf zet2dr var-11)
  455.             (setf zet2di var-12)
  456.             (setf asumdr var-13)
  457.             (setf asumdi var-14)
  458.             (setf bsumdr var-15)
  459.             (setf bsumdi var-16))
  460.          label210
  461.           (if (= kode 1) (go label220))
  462.           (setf str (+ zbr zet2dr))
  463.           (setf sti (+ zbi zet2di))
  464.           (setf rast (/ fn (zabs str sti)))
  465.           (setf str (* str rast rast))
  466.           (setf sti (* (- sti) rast rast))
  467.           (setf s1r (- str zet1dr))
  468.           (setf s1i (- sti zet1di))
  469.           (go label230)
  470.          label220
  471.           (setf s1r (- zet2dr zet1dr))
  472.           (setf s1i (- zet2di zet1di))
  473.          label230
  474.           (setf rs1 s1r)
  475.           (if (> (abs rs1) elim) (go label280))
  476.           (if (= kdflg 1) (setf iflag 2))
  477.           (if (< (abs rs1) alim) (go label240))
  478.           (setf aphi (zabs phidr phidi))
  479.           (setf aarg (zabs argdr argdi))
  480.           (setf rs1
  481.                   (- (+ rs1 (f2cl-lib:flog aphi))
  482.                      (* 0.25 (f2cl-lib:flog aarg))
  483.                      aic))
  484.           (if (> (abs rs1) elim) (go label280))
  485.           (if (= kdflg 1) (setf iflag 1))
  486.           (if (< rs1 0.0) (go label240))
  487.           (if (= kdflg 1) (setf iflag 3))
  488.          label240
  489.           (multiple-value-bind
  490.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  491.               (zairy argdr argdi 0 2 air aii nai idum)
  492.             (declare (ignore var-0 var-1 var-2 var-3))
  493.             (setf air var-4)
  494.             (setf aii var-5)
  495.             (setf nai var-6)
  496.             (setf idum var-7))
  497.           (multiple-value-bind
  498.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  499.               (zairy argdr argdi 1 2 dair daii ndai idum)
  500.             (declare (ignore var-0 var-1 var-2 var-3))
  501.             (setf dair var-4)
  502.             (setf daii var-5)
  503.             (setf ndai var-6)
  504.             (setf idum var-7))
  505.           (setf str (- (* dair bsumdr) (* daii bsumdi)))
  506.           (setf sti (+ (* dair bsumdi) (* daii bsumdr)))
  507.           (setf str (+ str (- (* air asumdr) (* aii asumdi))))
  508.           (setf sti (+ sti (+ (* air asumdi) (* aii asumdr))))
  509.           (setf ptr (- (* str phidr) (* sti phidi)))
  510.           (setf pti (+ (* str phidi) (* sti phidr)))
  511.           (setf s2r (- (* ptr csr) (* pti csi)))
  512.           (setf s2i (+ (* ptr csi) (* pti csr)))
  513.           (setf str (* (exp s1r) (f2cl-lib:fref cssr (iflag) ((1 3)))))
  514.           (setf s1r (* str (cos s1i)))
  515.           (setf s1i (* str (sin s1i)))
  516.           (setf str (- (* s2r s1r) (* s2i s1i)))
  517.           (setf s2i (+ (* s2r s1i) (* s2i s1r)))
  518.           (setf s2r str)
  519.           (if (/= iflag 1) (go label250))
  520.           (multiple-value-bind
  521.               (var-0 var-1 var-2 var-3 var-4)
  522.               (zuchk s2r s2i nw (f2cl-lib:fref bry (1) ((1 3))) tol)
  523.             (declare (ignore var-0 var-1 var-3 var-4))
  524.             (setf nw var-2))
  525.           (if (= nw 0) (go label250))
  526.           (setf s2r zeror)
  527.           (setf s2i zeroi)
  528.          label250
  529.           (if (<= yy 0.0) (setf s2i (- s2i)))
  530.           (f2cl-lib:fset (f2cl-lib:fref cyr (kdflg) ((1 2))) s2r)
  531.           (f2cl-lib:fset (f2cl-lib:fref cyi (kdflg) ((1 2))) s2i)
  532.           (setf c2r s2r)
  533.           (setf c2i s2i)
  534.           (setf s2r (* s2r (f2cl-lib:fref csrr (iflag) ((1 3)))))
  535.           (setf s2i (* s2i (f2cl-lib:fref csrr (iflag) ((1 3)))))
  536.           (setf s1r (f2cl-lib:fref yr (kk) ((1 n))))
  537.           (setf s1i (f2cl-lib:fref yi (kk) ((1 n))))
  538.           (if (= kode 1) (go label270))
  539.           (multiple-value-bind
  540.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  541.               (zs1s2 zrr zri s1r s1i s2r s2i nw asc alim iuf)
  542.             (declare (ignore var-0 var-1 var-7 var-8))
  543.             (setf s1r var-2)
  544.             (setf s1i var-3)
  545.             (setf s2r var-4)
  546.             (setf s2i var-5)
  547.             (setf nw var-6)
  548.             (setf iuf var-9))
  549.           (setf nz (f2cl-lib:int-add nz nw))
  550.          label270
  551.           (f2cl-lib:fset (f2cl-lib:fref yr (kk) ((1 n)))
  552.                          (+ (- (* s1r cspnr) (* s1i cspni)) s2r))
  553.           (f2cl-lib:fset (f2cl-lib:fref yi (kk) ((1 n)))
  554.                          (+ (* s1r cspni) (* s1i cspnr) s2i))
  555.           (setf kk (f2cl-lib:int-sub kk 1))
  556.           (setf cspnr (- cspnr))
  557.           (setf cspni (- cspni))
  558.           (setf str csi)
  559.           (setf csi (- csr))
  560.           (setf csr str)
  561.           (if (or (/= c2r 0.0) (/= c2i 0.0)) (go label255))
  562.           (setf kdflg 1)
  563.           (go label290)
  564.          label255
  565.           (if (= kdflg 2) (go label295))
  566.           (setf kdflg 2)
  567.           (go label290)
  568.          label280
  569.           (if (> rs1 0.0) (go label320))
  570.           (setf s2r zeror)
  571.           (setf s2i zeroi)
  572.           (go label250)
  573.          label290))
  574.       (setf k n)
  575.      label295
  576.       (setf il (f2cl-lib:int-sub n k))
  577.       (if (= il 0) (go end_label))
  578.       (setf s1r (f2cl-lib:fref cyr (1) ((1 2))))
  579.       (setf s1i (f2cl-lib:fref cyi (1) ((1 2))))
  580.       (setf s2r (f2cl-lib:fref cyr (2) ((1 2))))
  581.       (setf s2i (f2cl-lib:fref cyi (2) ((1 2))))
  582.       (setf csr (f2cl-lib:fref csrr (iflag) ((1 3))))
  583.       (setf ascle (f2cl-lib:fref bry (iflag) ((1 3))))
  584.       (setf fn
  585.               (coerce (the f2cl-lib:integer4 (f2cl-lib:int-add inu il))
  586.                       'double-float))
  587.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  588.                     ((> i il) nil)
  589.         (tagbody
  590.           (setf c2r s2r)
  591.           (setf c2i s2i)
  592.           (setf s2r (+ s1r (* (+ fn fnf) (- (* rzr c2r) (* rzi c2i)))))
  593.           (setf s2i (+ s1i (* (+ fn fnf) (+ (* rzr c2i) (* rzi c2r)))))
  594.           (setf s1r c2r)
  595.           (setf s1i c2i)
  596.           (setf fn (- fn 1.0))
  597.           (setf c2r (* s2r csr))
  598.           (setf c2i (* s2i csr))
  599.           (setf ckr c2r)
  600.           (setf cki c2i)
  601.           (setf c1r (f2cl-lib:fref yr (kk) ((1 n))))
  602.           (setf c1i (f2cl-lib:fref yi (kk) ((1 n))))
  603.           (if (= kode 1) (go label300))
  604.           (multiple-value-bind
  605.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  606.               (zs1s2 zrr zri c1r c1i c2r c2i nw asc alim iuf)
  607.             (declare (ignore var-0 var-1 var-7 var-8))
  608.             (setf c1r var-2)
  609.             (setf c1i var-3)
  610.             (setf c2r var-4)
  611.             (setf c2i var-5)
  612.             (setf nw var-6)
  613.             (setf iuf var-9))
  614.           (setf nz (f2cl-lib:int-add nz nw))
  615.          label300
  616.           (f2cl-lib:fset (f2cl-lib:fref yr (kk) ((1 n)))
  617.                          (+ (- (* c1r cspnr) (* c1i cspni)) c2r))
  618.           (f2cl-lib:fset (f2cl-lib:fref yi (kk) ((1 n)))
  619.                          (+ (* c1r cspni) (* c1i cspnr) c2i))
  620.           (setf kk (f2cl-lib:int-sub kk 1))
  621.           (setf cspnr (- cspnr))
  622.           (setf cspni (- cspni))
  623.           (if (>= iflag 3) (go label310))
  624.           (setf c2r (coerce (abs ckr) 'double-float))
  625.           (setf c2i (coerce (abs cki) 'double-float))
  626.           (setf c2m (max c2r c2i))
  627.           (if (<= c2m ascle) (go label310))
  628.           (setf iflag (f2cl-lib:int-add iflag 1))
  629.           (setf ascle (f2cl-lib:fref bry (iflag) ((1 3))))
  630.           (setf s1r (* s1r csr))
  631.           (setf s1i (* s1i csr))
  632.           (setf s2r ckr)
  633.           (setf s2i cki)
  634.           (setf s1r (* s1r (f2cl-lib:fref cssr (iflag) ((1 3)))))
  635.           (setf s1i (* s1i (f2cl-lib:fref cssr (iflag) ((1 3)))))
  636.           (setf s2r (* s2r (f2cl-lib:fref cssr (iflag) ((1 3)))))
  637.           (setf s2i (* s2i (f2cl-lib:fref cssr (iflag) ((1 3)))))
  638.           (setf csr (f2cl-lib:fref csrr (iflag) ((1 3))))
  639.          label310))
  640.       (go end_label)
  641.      label320
  642.       (setf nz -1)
  643.       (go end_label)
  644.      end_label
  645.       (return (values nil nil nil nil nil nil nil nil nz nil nil nil)))))
  646.  
  647.